#ifndef Kernel_MModel_h
#define Kernel_MModel_h

#include <set>

#include "Utility/Algebra/Config.h"
#include "Utility/Algebra/ODESystem.h"

#include "DataStructure/Common/Id.h"
#include "DataStructure/Common/Out.h"
#include "DataStructure/Common/Real.h"
#include "DataStructure/Input/Dynamic/InitConditionType.h"
#include "DataStructure/Input/Solver/AnaysisTypes.h"
#include "DataStructure/Input/Load/DLOADData.h"
#include "DataStructure/Result/Random/CROSS_PSDData.h"

#include "SGFem/Pre/Assembly/Assembly.h"    //为使用EigenSparse别名，临时从.cpp放置此处

#include "DBManager/DBServiceFactory.h"
#include "DBManager/Common/DomainRequest.h"

#include "Domain.h"
#include "KernelModule.h"
#include "SGFem/Task/Common/DynamicLoad.h"
#include "TimeIndex.h"

namespace Kernel
{
    /**
     * \addtogroup Kernel
     * @{
     */

    // Forward declarations
    class Variable;

    /**
     * @brief Represent the physical/mathematical model under consideration.
     *
     * MModel works as a facade to shield the upper algorithms, namely tasks, from the underlying
     * FEM components and data. It knows the discretized representations of the governing equations
     * and understands how to gather local contributions over nodes/elements to form the corresponding
     * algebra terms.
     *
     * Facade Pattern is used to provide a unified interface to a set of interfaces in a subsystem.
     * Facade defines a higher-level interface that makes the subsystem easier to use.
     *
     * @see GoF's Facade Pattern.
     *
     * For the parallel mode on the distributed processors, the whole mesh will be partitioned
     * into several domains. And then it will contains the domain associated with a particular process.
     *
     * In summary, MModel maintains the following responsibilities.
     * - Assemble matrix/vector from local contributions over the entities (nodes/elements).
     * - Provide access to the global fields, such as displacement, velocity and so on.
     * - Relay requirements, for example, postprocessing, to the underlying FEM components.
     *
     * @note For now, only structural analysis has been considered.
     *
     * @author nene
     * @date November, 2024.
     *
     */
    class KERNEL_EXPORT MModel : public SG::Algebra::SecondOrderODESystem
    {
      public:
        MModel (SG::DataStructure::FEM::SOLType analysisType, SG::DBManager::DBServiceFactorySP dbServiceFactory);
        virtual ~MModel ();

        using Real    = SG::DataStructure::Common::Real;
        using I64     = SG::DataStructure::Common::I64;
        using Id      = SG::DataStructure::Common::Id;
        using IdArray = SG::DataStructure::Common::IdArray;
        using RandomResultType = SG::DBManager::Post::DomainRequest::RandomResultType;


        /**
         * @brief    开辟矩阵和自由度空间
         *
         * 节点和自由度排序
         * 坐标系基矢量初始化等
         *
         * @param    subcaseIds 一组工况的集合， 这些工况应该具有相同的边界条件
         * @author   wanghuan
         * @date     2025.07.18
         */
        virtual void Initialize (const IdArray& subcaseIds);

        /** @name Assemble from local contributions over the entities (nodes/elements) */
        ///@{
        /**
         * @brief 组装 K M 以及静力载荷
         *
         * @param subcaseId        用于确定边界条件
         * @param loadSidSet       计算所有有效载荷数据
         */
        SG::DataStructure::Common::Status CalcMatrixEquation (int subcaseId);

        /**
         * @brief   计算有限元方程右端荷载向量   
         * 
         * @param    dbServiceFactory 
         * @param    loadId 
         * @param    f 
         * @author   Liu Xing
         * @date     2025.11.03
         */
        void ComputeLoadVector(const SG::DataStructure::Common::Id loadId, SG::FEMSolver::EigenVector& f) const;
        
        void ComputeStiffnessMatrix (SG::Algebra::SpMat_t& K);
        void ComputeMassMatrix (SG::Algebra::SpMat_t& M);
        void ComputeLumpedMassMatrix (SG::Algebra::Vec_t& M);

        /**
         * @brief    计算粘滞阻尼
         *
         * 瞬态响应分析中：C = B1 + B2 + G/W3*K + 1/W4*∑(GE*KE)；
         * 频率响应分析中：C = B1 + B2
         *
         * @param    C
         * @author   wjw
         * @date     2025.03.31
         */
        void                ComputeDampingMatrix (SG::Algebra::SpMat_t& C);
        SG::Algebra::Real_t ComputeDampingRatio (SG::Algebra::Real_t omega);

        /**
         * @brief    计算动力载荷的空间分布
         *
         * @param    subcaseId
         * @author   wanghuan
         * @date     2025.09.05
         */
        void ComputeFValue (const Id subcaseId);
        // void ComputeForce (SG::Algebra::Vec_t& F);
        ///@}

        /** @name Transient solver */
        ///@{
      private:
        void TimeDomainLoadInit (const SG::DataStructure::FEM::DLOADData& dloadData, SG::DataStructure::Common::Id loadSetId);
      public:
        SG::DataStructure::Common::Status QueryTimeStep (Id                                            tStepId,
                                                         std::vector<SG::DataStructure::Common::Real>& tCompute,
                                                         std::vector<SG::DataStructure::Common::Real>& tOut) const;

        void ComputeInitialCondition (Id                                       icId,
                                      SG::DataStructure::FEM::InitConditonType type,
                                      SG::Algebra::Vec_t&                      U,
                                      SG::Algebra::Vec_t&                      V) const;
        ///@}

        /** @name Frequency solver */
        ///@{
      private:
        void FreqDomainLoadInit (const SG::DataStructure::FEM::DLOADData& dloadData, SG::DataStructure::Common::Id loadSetId);
      public:
        /**
         * @brief    计算频响分析频率点
         *
         * @param    tbId      工况段指定的频率表id
         * @param    freSet
         * @return
         * @author   wanghuan
         * @date     2025.07.18
         */
        SG::DataStructure::Common::Status QueryFrenqencePoint (Id tbId, std::vector<SG::DataStructure::Common::Real>& _OUT freSet);

        /**
         * @brief    返回动刚度矩阵
         *
         * Kdyna = (1.0 + i G) * K +  i sum{Ge *Ke}
         * @param    Kdyna
         * @author   wanghuan
         * @date     2025.09.23
         */
        void ComputeDynaStiffness (SG::Algebra::ComplexSpMat_t& _OUT Kdyna);

        /// @brief 指定频率下的载荷向量
        /// @param t
        /// @param u
        /// @param v
        /// @param Pf
        void ComputeFreqForce (const SG::Algebra::Real_t& t,
                               const SG::Algebra::Vec_t&  u,
                               const SG::Algebra::Vec_t&  v,
                                   SG::Algebra::ComplexVec_t&        Pf);       
        ///@}

        void ComputeDerivatives (const SG::Algebra::Real_t& t,
                                 const SG::Algebra::Vec_t&  u,
                                 const SG::Algebra::Vec_t&  v,
                                 SG::Algebra::Vec_t&        a) override;
        void ComputeImplicitValue (const SG::Algebra::Real_t& t,
                                   const SG::Algebra::Vec_t&  u,
                                   const SG::Algebra::Vec_t&  v,
                                   const SG::Algebra::Vec_t&  a,
                                   SG::Algebra::Vec_t&        f) override;

        // 发现 频率和 时间域都使用了这个接口
        // 作为虚函数，会被时间积分器调用
        void ComputeExplicitValue (const SG::Algebra::Real_t& t,
                                   const SG::Algebra::Vec_t&  u,
                                   const SG::Algebra::Vec_t&  v,
                                   SG::Algebra::Vec_t&        g) override;

        void ComputeImplicitJacobian (const SG::Algebra::Real_t& t,
                                      const SG::Algebra::Vec_t&  u,
                                      const SG::Algebra::Vec_t&  v,
                                      const SG::Algebra::Vec_t&  a,
                                      SG::Algebra::SpMat_t&      dfdu,
                                      SG::Algebra::SpMat_t&      dfdv,
                                      SG::Algebra::SpMat_t&      dfda) override;
        // 对线性有限元此方法不需要干任何事情
        void ComputeExplicitJacobian (const SG::Algebra::Real_t& t,
                                      const SG::Algebra::Vec_t&  u,
                                      const SG::Algebra::Vec_t&  v,
                                      SG::Algebra::SpMat_t&      dgdu,
                                      SG::Algebra::SpMat_t&      dgdv) override {};


        /** @name Calculate post-analysis results */
        ///@{
        void ComputeOutputData (const TimeIndex& t);
        ///@}

        /** @name Field Manipulations */
        ///@{
        /// General methods to get/set the fields (Operate on the historical data).
        void GetNodalSolution (const Variable& variable, SG::Algebra::Vec_t& value, const TimeIndex& t = 0);
        void SetNodalSolution (const Variable& variable, SG::Algebra::Vec_t& value, const TimeIndex& t = 0);
        void GetNodalSolution (const Variable& variable, SG::Algebra::Mat_t& value, const TimeIndex& t = 0);
        void SetNodalSolution (const Variable& variable, SG::Algebra::Mat_t& value, const TimeIndex& t = 0);

        // Special methods dedicated to the 'displacement/velocity/acceleration/EigenVector/EigenValue' fields.
        void GetDisplacement (const Id subcaseId, const double t, _OUT SG::Algebra::Vec_t& V) const;
        void SetDisplacement (SG::Algebra::Vec_t V, const Id subcaseId, const double t);
        void SetDisplacement (const SG::Algebra::Vec_t& V, Id randomSid, SG::Algebra::Real_t t, RandomResultType random);
        void GetDisplacement (const Id subcaseId, const double t, _OUT SG::Algebra::ComplexVec_t& V) const;
        SG::DataStructure::Common::Status GetFreaRespMatrix (SG::DataStructure::Common::Real f, const IdArray& subcaseIds, SG::Algebra::ComplexMat_t& M);
        void SetDisplacement (SG::Algebra::ComplexVec_t& V, const Id subcaseId, const double t);

        void GetVelocity (_OUT SG::Algebra::Mat_t& V, const Id subcaseId, const TimeIndex& t = 0) const;
        void SetVelocity (SG::Algebra::Mat_t& V, const Id subcaseId, const TimeIndex& t = 0);
        void GetAcceleration (_OUT SG::Algebra::Mat_t& A, const Id subcaseId, const TimeIndex& t = 0) const;
        void SetAcceleration (SG::Algebra::Mat_t& A, const Id subcaseId, const TimeIndex& t = 0);
        ///@}

        void SetVelocity (SG::Algebra::Vec_t V, const Id subcaseId, const double t);
        void SetAcceleration (SG::Algebra::Vec_t V, const Id subcaseId, const double t);

        void GetVelocity (const Id subcaseId, const double t, _OUT SG::Algebra::Vec_t& V) const;
        void GetAcceleration (const Id subcaseId, const double t, _OUT SG::Algebra::Vec_t& V) const;

        /** @name Real eigenvector  compute*/
        ///@{
                /**
         * @brief    计算广义特征值，并保存特征值到数据库（因为保存特征值需要用到 K 和 M 的数据，因此在函数内部进行）
         * 
         * @param    ModesOrder 
         * @param    eigD
         * @param    eigV 
         * @author   wanghuan
         * @date     2025.10.21
         */
        void ComputeEigenValue (int ModesOrder, SG::Algebra::Vec_t& _OUT eigD, SG::Algebra::Mat_t& _OUT eigV);
        
        void SetEigenVector (SG::Algebra::Mat_t& eigenVectorM, SG::Algebra::Vec_t eigenValue, const Id subcaseId);
        void GetEigenVector (_OUT std::vector<Real>& V, const Id subcaseId, const TimeIndex& t = 0) const;

        void SetCrossPSD (const std::vector<SG::DataStructure::Post::CROSS_PSDData>& CROSSPSD, const Id randomSid, const SG::Algebra::Real_t freq);

        /**
         * @brief   : 模态分析中 保存某工况某阶特征值
         * @param   : subcaseId 工况ID
         * @param   : eigen KΦ=λMΦ里边的λ，就是特征值
         * @param   : mass  广义质量
         * @param   : stiff 广义刚度
         * @param   : t     阶
         * @param   : resFlag Residual vector augmentation flag
         * @param   : fldFlag Fluid modes flag
         * @author  : wxc
         * @date    : 2024.11.20
         */
        void SetEigenValue (const int t, const Real eigen, const Real mass, const Real stiff, const I64 resFlag = 0, const I64 fldFlag = 0);

        /**
         * @brief    Get the Eigen Value object
         *
         * @param    subcaseId  工况ID
         * @param    t 阶
         * @param    eigen
         * @param    omega
         * @param    freq
         * @param    mass
         * @param    stiff
         * @param    resFlag
         * @param    fldFlag
         * @author   wxc
         * @date     2024.11.28
         */
        void GetEigenValue (const I64  subcaseId,
                            const int  t,
                            _OUT Real& eigen,
                            _OUT Real& omega,
                            _OUT Real& freq,
                            _OUT Real& mass,
                            _OUT Real& stiff,
                            _OUT I64&  resFlag,
                            _OUT I64&  fldFlag) const;
        ///}@

        /** @name Complex eigenvector  compute*/
        ///@{
        void SetEigenVector (SG::Algebra::ComplexMat_t& complexMat, SG::Algebra::ComplexVec_t eigenValue, const Id subcaseId);
        void GetEigenVector (const Id subcaseId, const TimeIndex& t, _OUT std::vector<SG::Algebra::Complex_t>& V) const;

        void SetEigenValue (const int t, const SG::Algebra::Complex_t eigen, const Real freq, const Real damp, const Real spinSpeed);
        void GetEigenValue (const I64 subcaseId,
                            const int t,
                            _OUT SG::Algebra::Complex_t& eigen,
                            _OUT Real&                   freq,
                            _OUT Real&                   damp,
                            _OUT Real&                   spinSpeed) const;
        ///}@

        virtual void Finalize ();

        // 临时获取subcase， 未来删除
        SG::DBManager::DBServiceFactorySP GetDBServiceFactory () const
        {
            return m_dbServiceFactory;
        }

      protected:
        std::vector<std::shared_ptr<Domain>> myDomains;                     //< mesh partitions
        SG::DBManager::DBServiceFactorySP    m_dbServiceFactory = nullptr;  ///< Database service factory sharedPtr

        /**
         * 分析类型
         * 对应的Code参考：
         * https://nexus.hexagon.com/documentationcenter/zh-CH/bundle/MSC_Nastran_2024.1/page/Nastran_Combined_Book/dmap/datablocks/TOC.OFP.Tables.xhtml#id_d213908e417
         */
        SG::DataStructure::FEM::SOLType           m_analysisType = SG::DataStructure::FEM::SOLType::SOL_UNKNOWN;

        std::set<SG::DataStructure::Common::Id>   m_loadCases;
        /// @brief 时域有关载荷向量
        std::vector<TimeDynamicLoadP> m_timeDynaLoadSet;

        /// @brief 频域有关载荷向量
        std::vector<FreqDynamicLoadP> m_freqDynaLoadSet;
    };

    /** @}*/

}  // namespace Kernel

#endif  // Kernel_MModel_h
